! This module is an add on to BeamDyn to allow output of data at each blade node.
!
! Copyright (C) 2016-2017  Envision Energy USA, LTD   
!
MODULE BeamDyn_BldNdOuts_IO

   USE NWTC_Library
   USE NWTC_LAPACK
   USE BeamDyn_Subs
   USE BeamDyn_Types

   IMPLICIT NONE

   PRIVATE


   ! Outstanding issues
   !  1. Currently nothing is added to the summary file. If we add some output there, some changes either in the BeamDyn
   !     code (as distributed) will be needed, or changes here (reopen file and append).


   PUBLIC   :: BldNdOuts_InitOut
   PUBLIC   :: Calc_WriteBldNdOutput
   PUBLIC   :: BldNdOuts_SetParameters



      ! Parameters related to output length (number of characters allowed in the output data headers):

!   INTEGER(IntKi), PARAMETER      :: OutStrLenM1 = ChanLen-6    ! We are making these of the form B1Z###quantity, but note that the glue code adds the "B1" (turbine component) part

! ===================================================================================================
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary.
! ===================================================================================================
! This code was generated by Write_ChckOutLst.m at 12-Dec-2017 20:48:14.

     ! NOTES:
     !    (1) These parameters are in the order stored in "OutListParameters.xlsx"


     ! Sectional Loads:

   INTEGER(IntKi), PARAMETER      :: BldNd_FxL      =   1
   INTEGER(IntKi), PARAMETER      :: BldNd_FyL      =   2
   INTEGER(IntKi), PARAMETER      :: BldNd_FzL      =   3
   INTEGER(IntKi), PARAMETER      :: BldNd_MxL      =   4
   INTEGER(IntKi), PARAMETER      :: BldNd_MyL      =   5
   INTEGER(IntKi), PARAMETER      :: BldNd_MzL      =   6
   INTEGER(IntKi), PARAMETER      :: BldNd_Fxr      =   7
   INTEGER(IntKi), PARAMETER      :: BldNd_Fyr      =   8
   INTEGER(IntKi), PARAMETER      :: BldNd_Fzr      =   9
   INTEGER(IntKi), PARAMETER      :: BldNd_Mxr      =  10
   INTEGER(IntKi), PARAMETER      :: BldNd_Myr      =  11
   INTEGER(IntKi), PARAMETER      :: BldNd_Mzr      =  12


     ! Sectional Motions:

   INTEGER(IntKi), PARAMETER      :: BldNd_TDxr     =  13
   INTEGER(IntKi), PARAMETER      :: BldNd_TDyr     =  14
   INTEGER(IntKi), PARAMETER      :: BldNd_TDzr     =  15
   INTEGER(IntKi), PARAMETER      :: BldNd_RDxr     =  16
   INTEGER(IntKi), PARAMETER      :: BldNd_RDyr     =  17
   INTEGER(IntKi), PARAMETER      :: BldNd_RDzr     =  18
   INTEGER(IntKi), PARAMETER      :: BldNd_AbsXg    =  19
   INTEGER(IntKi), PARAMETER      :: BldNd_AbsYg    =  20
   INTEGER(IntKi), PARAMETER      :: BldNd_AbsZg    =  21
   INTEGER(IntKi), PARAMETER      :: BldNd_AbsXr    =  22
   INTEGER(IntKi), PARAMETER      :: BldNd_AbsYr    =  23
   INTEGER(IntKi), PARAMETER      :: BldNd_AbsZr    =  24
   INTEGER(IntKi), PARAMETER      :: BldNd_TVxg     =  25
   INTEGER(IntKi), PARAMETER      :: BldNd_TVyg     =  26
   INTEGER(IntKi), PARAMETER      :: BldNd_TVzg     =  27
   INTEGER(IntKi), PARAMETER      :: BldNd_TVxl     =  28
   INTEGER(IntKi), PARAMETER      :: BldNd_TVyl     =  29
   INTEGER(IntKi), PARAMETER      :: BldNd_TVzl     =  30
   INTEGER(IntKi), PARAMETER      :: BldNd_TVxr     =  31
   INTEGER(IntKi), PARAMETER      :: BldNd_TVyr     =  32
   INTEGER(IntKi), PARAMETER      :: BldNd_TVzr     =  33
   INTEGER(IntKi), PARAMETER      :: BldNd_RVxg     =  34
   INTEGER(IntKi), PARAMETER      :: BldNd_RVyg     =  35
   INTEGER(IntKi), PARAMETER      :: BldNd_RVzg     =  36
   INTEGER(IntKi), PARAMETER      :: BldNd_RVxl     =  37
   INTEGER(IntKi), PARAMETER      :: BldNd_RVyl     =  38
   INTEGER(IntKi), PARAMETER      :: BldNd_RVzl     =  39
   INTEGER(IntKi), PARAMETER      :: BldNd_RVxr     =  40
   INTEGER(IntKi), PARAMETER      :: BldNd_RVyr     =  41
   INTEGER(IntKi), PARAMETER      :: BldNd_RVzr     =  42
   INTEGER(IntKi), PARAMETER      :: BldNd_TAxl     =  43
   INTEGER(IntKi), PARAMETER      :: BldNd_TAyl     =  44
   INTEGER(IntKi), PARAMETER      :: BldNd_TAzl     =  45
   INTEGER(IntKi), PARAMETER      :: BldNd_TAxr     =  46
   INTEGER(IntKi), PARAMETER      :: BldNd_TAyr     =  47
   INTEGER(IntKi), PARAMETER      :: BldNd_TAzr     =  48
   INTEGER(IntKi), PARAMETER      :: BldNd_RAxl     =  49
   INTEGER(IntKi), PARAMETER      :: BldNd_RAyl     =  50
   INTEGER(IntKi), PARAMETER      :: BldNd_RAzl     =  51
   INTEGER(IntKi), PARAMETER      :: BldNd_RAxr     =  52
   INTEGER(IntKi), PARAMETER      :: BldNd_RAyr     =  53
   INTEGER(IntKi), PARAMETER      :: BldNd_RAzr     =  54

     ! Applied Loads:

   INTEGER(IntKi), PARAMETER      :: BldNd_PFxL     =  55
   INTEGER(IntKi), PARAMETER      :: BldNd_PFyL     =  56
   INTEGER(IntKi), PARAMETER      :: BldNd_PFzL     =  57
   INTEGER(IntKi), PARAMETER      :: BldNd_PMxL     =  58
   INTEGER(IntKi), PARAMETER      :: BldNd_PMyL     =  59
   INTEGER(IntKi), PARAMETER      :: BldNd_PMzL     =  60
   INTEGER(IntKi), PARAMETER      :: BldNd_DFxL     =  61
   INTEGER(IntKi), PARAMETER      :: BldNd_DFyL     =  62
   INTEGER(IntKi), PARAMETER      :: BldNd_DFzL     =  63
   INTEGER(IntKi), PARAMETER      :: BldNd_DMxL     =  64
   INTEGER(IntKi), PARAMETER      :: BldNd_DMyL     =  65
   INTEGER(IntKi), PARAMETER      :: BldNd_DMzL     =  66
   INTEGER(IntKi), PARAMETER      :: BldNd_DFxR     =  67
   INTEGER(IntKi), PARAMETER      :: BldNd_DFyR     =  68
   INTEGER(IntKi), PARAMETER      :: BldNd_DFzR     =  69
   INTEGER(IntKi), PARAMETER      :: BldNd_DMxR     =  70
   INTEGER(IntKi), PARAMETER      :: BldNd_DMyR     =  71
   INTEGER(IntKi), PARAMETER      :: BldNd_DMzR     =  72


     ! Sectional Partial Loads (debugging):

   INTEGER(IntKi), PARAMETER      :: BldNd_FFbxl    =  73
   INTEGER(IntKi), PARAMETER      :: BldNd_FFbyl    =  74
   INTEGER(IntKi), PARAMETER      :: BldNd_FFbzl    =  75
   INTEGER(IntKi), PARAMETER      :: BldNd_FFbxr    =  76
   INTEGER(IntKi), PARAMETER      :: BldNd_FFbyr    =  77
   INTEGER(IntKi), PARAMETER      :: BldNd_FFbzr    =  78
   INTEGER(IntKi), PARAMETER      :: BldNd_MFbxl    =  79
   INTEGER(IntKi), PARAMETER      :: BldNd_MFbyl    =  80
   INTEGER(IntKi), PARAMETER      :: BldNd_MFbzl    =  81
   INTEGER(IntKi), PARAMETER      :: BldNd_MFbxr    =  82
   INTEGER(IntKi), PARAMETER      :: BldNd_MFbyr    =  83
   INTEGER(IntKi), PARAMETER      :: BldNd_MFbzr    =  84
   INTEGER(IntKi), PARAMETER      :: BldNd_FFcxl    =  85
   INTEGER(IntKi), PARAMETER      :: BldNd_FFcyl    =  86
   INTEGER(IntKi), PARAMETER      :: BldNd_FFczl    =  87
   INTEGER(IntKi), PARAMETER      :: BldNd_FFcxr    =  88
   INTEGER(IntKi), PARAMETER      :: BldNd_FFcyr    =  89
   INTEGER(IntKi), PARAMETER      :: BldNd_FFczr    =  90
   INTEGER(IntKi), PARAMETER      :: BldNd_MFcxl    =  91
   INTEGER(IntKi), PARAMETER      :: BldNd_MFcyl    =  92
   INTEGER(IntKi), PARAMETER      :: BldNd_MFczl    =  93
   INTEGER(IntKi), PARAMETER      :: BldNd_MFcxr    =  94
   INTEGER(IntKi), PARAMETER      :: BldNd_MFcyr    =  95
   INTEGER(IntKi), PARAMETER      :: BldNd_MFczr    =  96
   INTEGER(IntKi), PARAMETER      :: BldNd_FFdxl    =  97
   INTEGER(IntKi), PARAMETER      :: BldNd_FFdyl    =  98
   INTEGER(IntKi), PARAMETER      :: BldNd_FFdzl    =  99
   INTEGER(IntKi), PARAMETER      :: BldNd_FFdxr    = 100
   INTEGER(IntKi), PARAMETER      :: BldNd_FFdyr    = 101
   INTEGER(IntKi), PARAMETER      :: BldNd_FFdzr    = 102
   INTEGER(IntKi), PARAMETER      :: BldNd_MFdxl    = 103
   INTEGER(IntKi), PARAMETER      :: BldNd_MFdyl    = 104
   INTEGER(IntKi), PARAMETER      :: BldNd_MFdzl    = 105
   INTEGER(IntKi), PARAMETER      :: BldNd_MFdxr    = 106
   INTEGER(IntKi), PARAMETER      :: BldNd_MFdyr    = 107
   INTEGER(IntKi), PARAMETER      :: BldNd_MFdzr    = 108
   INTEGER(IntKi), PARAMETER      :: BldNd_FFgxl    = 109
   INTEGER(IntKi), PARAMETER      :: BldNd_FFgyl    = 110
   INTEGER(IntKi), PARAMETER      :: BldNd_FFgzl    = 111
   INTEGER(IntKi), PARAMETER      :: BldNd_FFgxr    = 112
   INTEGER(IntKi), PARAMETER      :: BldNd_FFgyr    = 113
   INTEGER(IntKi), PARAMETER      :: BldNd_FFgzr    = 114
   INTEGER(IntKi), PARAMETER      :: BldNd_MFgxl    = 115
   INTEGER(IntKi), PARAMETER      :: BldNd_MFgyl    = 116
   INTEGER(IntKi), PARAMETER      :: BldNd_MFgzl    = 117
   INTEGER(IntKi), PARAMETER      :: BldNd_MFgxr    = 118
   INTEGER(IntKi), PARAMETER      :: BldNd_MFgyr    = 119
   INTEGER(IntKi), PARAMETER      :: BldNd_MFgzr    = 120
   INTEGER(IntKi), PARAMETER      :: BldNd_FFixl    = 121
   INTEGER(IntKi), PARAMETER      :: BldNd_FFiyl    = 122
   INTEGER(IntKi), PARAMETER      :: BldNd_FFizl    = 123
   INTEGER(IntKi), PARAMETER      :: BldNd_FFixr    = 124
   INTEGER(IntKi), PARAMETER      :: BldNd_FFiyr    = 125
   INTEGER(IntKi), PARAMETER      :: BldNd_FFizr    = 126
   INTEGER(IntKi), PARAMETER      :: BldNd_MFixl    = 127
   INTEGER(IntKi), PARAMETER      :: BldNd_MFiyl    = 128
   INTEGER(IntKi), PARAMETER      :: BldNd_MFizl    = 129
   INTEGER(IntKi), PARAMETER      :: BldNd_MFixr    = 130
   INTEGER(IntKi), PARAMETER      :: BldNd_MFiyr    = 131
   INTEGER(IntKi), PARAMETER      :: BldNd_MFizr    = 132

   ! The maximum number of output channels which can be output by the code.
   INTEGER(IntKi), PARAMETER, PUBLIC      :: BldNd_MaxOutPts      = 132

!End of code generated by Matlab script
! ===================================================================================================

CONTAINS
!----------------------------------------------------------------------------------------------------------------------------------

!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine populates the headers with the blade node outputs.  The iteration cycle is blade:node:channel (channel iterated
!! fastest).  If this iteration order is changed, it should be changed in the Calc_WriteBldNdOutput routine as well.
SUBROUTINE BldNdOuts_InitOut( InitOut, p, ErrStat, ErrMsg )


   TYPE(BD_InitOutputType),      INTENT(INOUT)  :: InitOut                          ! output data
   TYPE(BD_ParameterType),       INTENT(IN   )  :: p                                ! The module parameters
   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat                          ! The error status code
   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg                           ! The error message, if an error occurred

   INTEGER(IntKi)                               :: ErrStat2                         ! Error status code
   INTEGER(IntKi)                               :: INDX                             ! Index count within WriteOutput
   INTEGER(IntKi)                               :: IdxNode                          ! Counter to the blade node we ae on
   INTEGER(IntKi)                               :: IdxChan                          ! Counter to the channel we are outputting.
   CHARACTER(1)                                 :: ChanPrefix                       ! Name prefix (B#_ -- note that the B# part is added in FAST, not here)
   CHARACTER(4), ALLOCATABLE                    :: DistStr(:)                       ! Array of prefix (Z######y)
   CHARACTER(3)                                 :: TmpChar                          ! Temporary char array to hold the node digits (3 places only!!!!)
   CHARACTER(*), PARAMETER                      :: RoutineName = ('BldNdOuts_InitOut')


         ! Initialize some things
      ErrMsg = ''
      ErrStat = ErrID_None


      IF ( p%BldNd_NumOuts == 0 ) THEN
         return
      ENDIF

         ! create the channel names using the z-coordinate of the beam in mm
      ALLOCATE( DistStr(size(p%BldNd_BlOutNd)), STAT=ErrStat2)
         IF (ErrStat2 /= 0) THEN
            call SetErrStat(ErrID_Fatal, 'Error allocating DistStr array.', ErrStat, ErrMsg, RoutineName)
            return
         END IF

         ! Warn if we will run into issues with more than 999 nodes.
      IF (p%node_total > 999 ) CALL SetErrStat(ErrID_Severe,'More than 999 blade nodes in use.  Output channel headers will not '// &
            'correctly reflect blade stations beyond 999. Modifications to the variable ChanLen in FAST are required.',ErrStat,ErrMsg,RoutineName)

      DO IdxNode=1,size(p%BldNd_BlOutNd)
            ! Create the name prefix:
         WRITE (TmpChar,'(I3.3)')  IdxNode
         DistStr(IdxNode) = 'N' // TmpChar
      END DO


         ! Populate the header an unit lines for all blades and nodes
         ! First set a counter so we know where in the output array we are in
         ! NOTE: we populate invalid names as well (some names are not valid outputs for certain configurations).  That means we will have zeros in those values.
      INDX = p%NumOuts + 1       ! p%NumOuts is the number of outputs from the normal BeamDyn output.  The WriteOutput array is sized to p%NumOuts + num(BldNdOuts)


!      ChanPrefix = '_' !newer names have underscore character to deliniate between sections
      ChanPrefix = ''
      ! Loop through all the outputs we requested here:
      DO IdxChan=1,p%BldNd_NumOuts

         DO IdxNode=1,size(p%BldNd_BlOutNd)

               ! Now write to the header
            InitOut%WriteOutputHdr(INDX) = TRIM(ChanPrefix) // TRIM(DistStr(IdxNode)) // p%BldNd_OutParam(IdxChan)%Name
            InitOut%WriteOutputUnt(INDX) = p%BldNd_OutParam(IdxChan)%Units

               ! Increment the index to the Header arrays
            INDX = INDX + 1

         ENDDO

      ENDDO

      IF (ALLOCATED(DistStr)) DEALLOCATE(DistStr)


END SUBROUTINE BldNdOuts_InitOut

!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine populates the headers with the blade node outputs.  The iteration cycle is blade:node:channel (channel iterated
!! fastest).  If this iteration order is changed, it should be changed in the Calc_WriteBldNdOutput routine as well.
SUBROUTINE Calc_WriteBldNdOutput( p, OtherState, m, y, ErrStat, ErrMsg )
   TYPE(BD_ParameterType),       INTENT(IN   )  :: p                                ! The module parameters
   TYPE(BD_OtherStateType),      INTENT(IN   )  :: OtherState                       ! Other states at time t
   TYPE(BD_MiscVarType),         INTENT(INOUT)  :: m                                ! misc variables
   TYPE(BD_OutputType),          INTENT(INOUT)  :: y                                ! outputs
   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat                          ! The error status code
   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg                           ! The error message, if an error occurred

      ! local variables
   INTEGER(IntKi)                               :: IdxOutList                       ! Index within WriteOutput
   INTEGER(IntKi)                               :: IdxChan                          ! Counter to the channel we are outputting.
   CHARACTER(*), PARAMETER                      :: RoutineName = 'Calc_WriteBldNdOutput'


      ! temporary variables for calculation
   INTEGER(IntKi)                               :: idx_node_in_elem                 ! node on current element
   INTEGER(IntKi)                               :: nelem                            ! current element
   INTEGER(IntKi)                               :: idx_node                         ! Counter to the blade node we are on
   INTEGER(IntKi)                               :: compIndx                         ! index for array component (x,y,z)
   REAL(BDKi)                                   :: Tmp33a(3,3)                      ! Temporary 3x4 for orientation calcs
   REAL(BDKi)                                   :: Tmp33b(3,3)                      ! Temporary 3x4 for orientation calcs

   REAL(BDKi)                                   :: WM_ParamRD(3)                    ! Wiener Milenkovic parameters for current node, in Global coordinates
   REAL(BDKi)                                   :: temp_vec(3)                      ! temporary vector for orientation info.
   REAL(BDKi)                                   :: temp_vec2(3)                     ! temporary vector for orientation info.
   REAL(BDKi)                                   :: temp_vec3(3)                     ! temporary vector for orientation info.
   REAL(BDKi)                                   :: d_ref(3)                         ! root displacement
   REAL(BDKi)                                   :: d(3)                             ! displacement

      ! WM param finding
   REAL(BDKi)                                   :: RootRelOrient(3,3)

         ! Error handling
   INTEGER(IntKi)          :: ErrStat2                     ! Temporary Error status
   CHARACTER(ErrMsgLen)    :: ErrMsg2                      ! Temporary Error message 


      ! Initialize some things
   ErrMsg = ''
   ErrStat = ErrID_None

   IF ( p%BldNd_NumOuts == 0 ) THEN
      return
   ENDIF


      ! Set the root rotation DCM relative to the reference.
      ! NOTE: the orientations used in this routine are DCM's.  These are directly from the mesh.
   call LAPACK_GEMM('T', 'N', 1.0_BDKi, m%u2%RootMotion%Orientation(:,:,1), m%u2%RootMotion%RefOrientation(:,:,1), 0.0_BDKi, RootRelOrient,   ErrStat2, ErrMsg2 )


      ! Loop over the channel sets
   DO IdxChan=1,p%BldNd_NumOuts


            ! Case to assign output to this channel and populate based on Indx value (this indicates what the channel is)
            ! Logic and mathematics used here come from Calc_WriteOutput

         SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output

            CASE (0) ! This is an invalid channel so we'll just skip it
               CYCLE
               
            !----------------------------------------
            ! Sectional translational locations and deflections (relative to the undeflected position) expressed in g
            CASE (BldNd_AbsXg,BldNd_AbsYg,BldNd_AbsZg)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_AbsXg)
                  compIndx = 1
               CASE (BldNd_AbsYg)
                  compIndx = 2
               CASE (BldNd_AbsZg)
                  compIndx = 3
               END SELECT
               
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  d     = y%BldMotion%TranslationDisp(:, idx_node)
                  d_ref = y%BldMotion%Position(       :, idx_node)
                     ! For actual global location
                  temp_vec = d + d_ref
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
               ENDDO
            !----------------------------------------
            ! Sectional translational locations and deflections (relative to the undeflected position) expressed in r
            CASE (BldNd_TDxr,BldNd_TDyr,BldNd_TDzr,BldNd_AbsXr,BldNd_AbsYr,BldNd_AbsZr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  d     = y%BldMotion%TranslationDisp(:, idx_node) - m%u2%RootMotion%TranslationDisp(:,1)
                  d_ref = y%BldMotion%Position(       :, idx_node) - m%u2%RootMotion%Position(       :,1)
                     ! For relative change in location
                  temp_vec2 = d + d_ref - matmul( RootRelOrient, d_ref )
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2)
                     ! For actual location relative to root
                  temp_vec2 = d + d_ref
                  temp_vec3 = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2)
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_TDxr)
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_TDyr)
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_TDzr)
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                     CASE (BldNd_AbsXr)
                        y%WriteOutput( IdxOutList ) = temp_vec3(1)
                     CASE (BldNd_AbsYr)
                        y%WriteOutput( IdxOutList ) = temp_vec3(2)
                     CASE (BldNd_AbsZr)
                        y%WriteOutput( IdxOutList ) = temp_vec3(3)
                  END SELECT
               ENDDO



            !----------------------------------------
            ! Rotational displacements as W-M parameters
            CASE ( BldNd_RDxr, BldNd_RDyr, BldNd_RDzr )
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_RDxr)
                  compIndx = 1
               CASE (BldNd_RDyr)
                  compIndx = 2
               CASE (BldNd_RDzr)
                  compIndx = 3
               END SELECT

               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output

                  !-------------------------
!FIXME: we are not trapping errors here.  Do we need to?
                  ! Sectional angular/rotational deflection Wiener-Milenkovic parameter (relative to the undeflected orientation) expressed in r
                  call LAPACK_GEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,idx_node), RootRelOrient,   0.0_BDKi, Tmp33b, ErrStat2, ErrMsg2 )
                  call LAPACK_GEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation(   :,:,idx_node), Tmp33b,          0.0_BDKi, Tmp33a,   ErrStat2, ErrMsg2 )
                  call BD_CrvExtractCrv(Tmp33a,temp_vec2, ErrStat2, ErrMsg2) ! temp_vec2 = the Wiener-Milenkovic parameters of the node's angular/rotational defelctions
                  WM_ParamRD = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2)    ! Rotate the parameters to the correct coordinate system for output

                  y%WriteOutput( IdxOutList ) = WM_ParamRD(compIndx)
               END DO


            !----------------------------------------
            ! Translational Velocities, global frame
            CASE (BldNd_TVxg)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  y%WriteOutput( IdxOutList ) = y%BldMotion%TranslationVel(1,idx_node) 
               ENDDO
            CASE (BldNd_TVyg)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  y%WriteOutput( IdxOutList ) = y%BldMotion%TranslationVel(2,idx_node) 
               ENDDO
            CASE (BldNd_TVzg)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  y%WriteOutput( IdxOutList ) = y%BldMotion%TranslationVel(3,idx_node) 
               ENDDO
            !----------------------------------------
            ! Rotational Velocities, global frame
            CASE (BldNd_RVxg)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  y%WriteOutput( IdxOutList ) = y%BldMotion%RotationVel(1,idx_node) * R2D
               ENDDO
            CASE (BldNd_RVyg)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  y%WriteOutput( IdxOutList ) = y%BldMotion%RotationVel(2,idx_node) * R2D
               ENDDO
            CASE (BldNd_RVzg)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  y%WriteOutput( IdxOutList ) = y%BldMotion%RotationVel(3,idx_node) * R2D
               ENDDO


            !----------------------------------------
            ! Translational Velocities, local frame
            CASE (BldNd_TVxl,BldNd_TVyl,BldNd_TVzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_TVxl)
                  compIndx = 1
               CASE (BldNd_TVyl)
                  compIndx = 2
               CASE (BldNd_TVzl)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node),y%BldMotion%TranslationVel(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
               ENDDO
            !----------------------------------------
            ! Rotational Velocities, local frame
            CASE (BldNd_RVxl,BldNd_RVyl,BldNd_RVzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_RVxl)
                  compIndx = 1
               CASE (BldNd_RVyl)
                  compIndx = 2
               CASE (BldNd_RVzl)
                  compIndx = 3
               END SELECT
               
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node),y%BldMotion%RotationVel(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)*R2D
               ENDDO


            !----------------------------------------
            ! Translational Velocities, root frame
            CASE (BldNd_TVxr,BldNd_TVyr,BldNd_TVzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_TVxr)
                  compIndx = 1
               CASE (BldNd_TVyr)
                  compIndx = 2
               CASE (BldNd_TVzr)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),y%BldMotion%TranslationVel(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
               ENDDO
            !----------------------------------------
            ! Rotational Velocities, root frame
            CASE (BldNd_RVxr,BldNd_RVyr,BldNd_RVzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_RVxr)
                  compIndx = 1
               CASE (BldNd_RVyr)
                  compIndx = 2
               CASE (BldNd_RVzr)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),y%BldMotion%RotationVel(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)*R2D
               ENDDO


            !----------------------------------------
            ! Translational Accelerations, local frame
            CASE (BldNd_TAxl, BldNd_TAyl, BldNd_TAzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_TAxl)
                  compIndx = 1
               CASE (BldNd_TAyl)
                  compIndx = 2
               CASE (BldNd_TAzl)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node),y%BldMotion%TranslationAcc(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
               ENDDO
            !----------------------------------------
            ! Rotational Accelerations, local frame
            CASE (BldNd_RAxl, BldNd_RAyl, BldNd_RAzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_RAxl)
                  compIndx = 1
               CASE (BldNd_RAyl)
                  compIndx = 2
               CASE (BldNd_RAzl)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node),y%BldMotion%RotationAcc(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)*R2D
               ENDDO


            !----------------------------------------
            ! Translational Accelerations, root frame
            CASE (BldNd_TAxr, BldNd_TAyr, BldNd_TAzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_TAxr)
                  compIndx = 1
               CASE (BldNd_TAyr)
                  compIndx = 2
               CASE (BldNd_TAzr)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),y%BldMotion%TranslationAcc(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
               ENDDO
            !----------------------------------------
            ! Rotational Accelerations, root frame
            CASE (BldNd_RAxr, BldNd_RAyr, BldNd_RAzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_RAxr)
                  compIndx = 1
               CASE (BldNd_RAyr)
                  compIndx = 2
               CASE (BldNd_RAzr)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes   ! Index to current output
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),y%BldMotion%RotationAcc(:,idx_node))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)*R2D
               ENDDO

            !----------------------------------------
            ! Applied point forces, local coordinate system. Not used when coupled to FAST.
            CASE (BldNd_PFxl,BldNd_PFyl,BldNd_PFzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_PFxl)
                  compIndx = 1
               CASE (BldNd_PFyl)
                  compIndx = 2
               CASE (BldNd_PFzl)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_FE) THEN
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%u2%PointLoad%Force( :,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else ! we need to do a mesh mapping first
!FIXME: this is not implemented yet.  
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     y%WriteOutput( IdxOutList ) = 0.0_ReKi
                  ENDDO
               end if
               
            !----------------------------------------
            ! Applied point moments, local coordinate system. Not used when coupled to FAST.
            CASE (BldNd_PMxl, BldNd_PMyl, BldNd_PMzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_PMxl)
                  compIndx = 1
               CASE (BldNd_PMyl)
                  compIndx = 2
               CASE (BldNd_PMzl)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                  y%WriteOutput( IdxOutList ) = 0.0_ReKi
               ENDDO
               if (p%BldMotionNodeLoc == BD_MESH_FE) THEN
                  DO idx_node=1,y%BldMotion%NNodes
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%u2%PointLoad%Moment( :,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else ! we need to do a mesh mapping first
!FIXME: this is not implemented yet.  
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     y%WriteOutput( IdxOutList ) = 0.0_ReKi
                  ENDDO
               end if

            !----------------------------------------
            ! Applied distributed forces (from AD15 when coupled to FAST), local frame
            CASE (BldNd_DFxl,BldNd_DFyl,BldNd_DFzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_DFxl)
                  compIndx = 1
               CASE (BldNd_DFyl)
                  compIndx = 2
               CASE (BldNd_DFzl)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_QP) THEN ! If we are on the quadrature points, the input and output meshes are siblings
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%u2%DistrLoad%Force( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%u_DistrLoad_at_y%Force( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if
            !----------------------------------------
            ! Applied distributed moments (from AD15 when coupled to FAST), local frame
            CASE (BldNd_DMxl,BldNd_DMyl,BldNd_DMzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_DMxl)
                  compIndx = 1
               CASE (BldNd_DMyl)
                  compIndx = 2
               CASE (BldNd_DMzl)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_QP) THEN ! If we are on the quadrature points, the input and output meshes are siblings
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%u2%DistrLoad%Moment( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%u_DistrLoad_at_y%Moment( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if

            !----------------------------------------
            ! Applied distributed forces (from AD15 when coupled to FAST), root frame
            CASE (BldNd_DFxr,BldNd_DFyr,BldNd_DFzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_DFxr)
                  compIndx = 1
               CASE (BldNd_DFyr)
                  compIndx = 2
               CASE (BldNd_DFzr)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_QP) THEN ! If we are on the quadrature points, the input and output meshes are siblings
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%u2%DistrLoad%Force( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%u_DistrLoad_at_y%Force( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if
            !----------------------------------------
            ! Applied distributed moments (from AD15 when coupled to FAST), root frame
            CASE (BldNd_DMxr,BldNd_DMyr,BldNd_DMzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_DMxr)
                  compIndx = 1
               CASE (BldNd_DMyr)
                  compIndx = 2
               CASE (BldNd_DMzr)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_QP) THEN ! If we are on the quadrature points, the input and output meshes are siblings
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%u2%DistrLoad%Moment( :,idx_node))
                        y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%u_DistrLoad_at_y%Moment( :,idx_node))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if


            !----------------------------------------
            ! Internal Forces, local
            CASE (BldNd_Fxl,BldNd_Fyl,BldNd_Fzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_Fxl)
                  compIndx = 1
               CASE (BldNd_Fyl)
                  compIndx = 2
               CASE (BldNd_Fzl)
                  compIndx = 3
               END SELECT
               
               if (p%BldMotionNodeLoc == BD_MESH_QP) then
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%BldInternalForceQP(1:3,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%BldInternalForceFE(1:3,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if
               
            !----------------------------------------
            ! Internal Moments, local
            CASE (BldNd_Mxl,BldNd_Myl,BldNd_Mzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_Mxl)
                  compIndx = 1
               CASE (BldNd_Myl)
                  compIndx = 2
               CASE (BldNd_Mzl)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_QP) then
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%BldInternalForceQP(4:6,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), m%BldInternalForceFE(4:6,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if

            !----------------------------------------
            ! Internal Forces, root frame
            CASE (BldNd_Fxr,BldNd_Fyr,BldNd_Fzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_Fxr)
                  compIndx = 1
               CASE (BldNd_Fyr)
                  compIndx = 2
               CASE (BldNd_Fzr)
                  compIndx = 3
               END SELECT
               if (p%BldMotionNodeLoc == BD_MESH_QP) then
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%BldInternalForceQP(1:3,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%BldInternalForceFE(1:3,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if

            !----------------------------------------
            ! Internal Moments, root frame
            CASE (BldNd_Mxr,BldNd_Myr,BldNd_Mzr)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )
               CASE (BldNd_Mxr)
                  compIndx = 1
               CASE (BldNd_Myr)
                  compIndx = 2
               CASE (BldNd_Mzr)
                  compIndx = 3
               END SELECT
               
               if (p%BldMotionNodeLoc == BD_MESH_QP) then
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%BldInternalForceQP(4:6,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               else
                  DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                     IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                     temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%BldInternalForceFE(4:6,p%NdIndxInverse(idx_node)))
                     y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
                  ENDDO
               end if
            
!>>> The remaining loads outputs are for debugging, and are not valid with BD_MESH_FE. We cannot get here in that case.
            !----------------------------------------
            ! Internal forces from CalcOutput, local frame
            CASE (BldNd_FFbxl,BldNd_FFbyl,BldNd_FFbzl)
               SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx ) 
               CASE (BldNd_FFbxl)
                  compIndx = 1
               CASE (BldNd_FFbyl)
                  compIndx = 2
               CASE (BldNd_FFbzl)
                  compIndx = 3
               END SELECT
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fb(1:3,idx_node_in_elem,nelem)))
                  y%WriteOutput( IdxOutList ) = temp_vec(compIndx)
               ENDDO
            CASE (BldNd_MFbxl,BldNd_MFbyl,BldNd_MFbzl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fb(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFbxl)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFbyl)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFbzl)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            CASE (BldNd_FFcxl,BldNd_FFcyl,BldNd_FFczl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fc(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFcxl)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFcyl)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFczl)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
              ENDDO
            CASE (BldNd_MFcxl,BldNd_MFcyl,BldNd_MFczl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fc(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFcxl)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFcyl)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFczl)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            CASE (BldNd_FFdxl,BldNd_FFdyl,BldNd_FFdzl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fd(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFdxl)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFdyl)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFdzl)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFdxl,BldNd_MFdyl,BldNd_MFdzl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fd(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFdxl)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFdyl)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFdzl)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            CASE (BldNd_FFgxl,BldNd_FFgyl,BldNd_FFgzl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fg(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFgxl)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFgyl)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFgzl)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFgxl,BldNd_MFgyl,BldNd_MFgzl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fg(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFgxl)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFgyl)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFgzl)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO



            !----------------------------------------
            ! Internal forces from CalcOutput, local frame
            CASE (BldNd_FFbxr,BldNd_FFbyr,BldNd_FFbzr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fb(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFbxr)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFbyr)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFbzr)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFbxr,BldNd_MFbyr,BldNd_MFbzr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fb(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFbxr)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFbyr)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFbzr)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            CASE (BldNd_FFcxr,BldNd_FFcyr,BldNd_FFczr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fc(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFcxr)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFcyr)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFczr)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFcxr,BldNd_MFcyr,BldNd_MFczr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fc(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFcxr)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFcyr)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFczr)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            CASE (BldNd_FFdxr,BldNd_FFdyr,BldNd_FFdzr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fd(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFdxr)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFdyr)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFdzr)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFdxr,BldNd_MFdyr,BldNd_MFdzr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fd(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFdxr)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFdyr)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFdzr)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            CASE (BldNd_FFgxr,BldNd_FFgyr,BldNd_FFgzr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fg(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFgxr)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFgyr)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFgzr)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFgxr,BldNd_MFgyr,BldNd_MFgzr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), m%qp%Fg(4:6,idx_node_in_elem,nelem))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFgxr)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFgyr)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFgzr)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO


            !----------------------------------------
            ! Inertial force from UpdateStates (Includes a few other terms), local frame
            CASE (BldNd_FFixl,BldNd_FFiyl,BldNd_FFizl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fi(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFixl)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFiyl)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFizl)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFixl,BldNd_MFiyl,BldNd_MFizl)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(y%BldMotion%Orientation(:,:,idx_node), MATMUL(OtherState%GlbRot,m%qp%Fi(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFixl)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFiyl)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFizl)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO

            !----------------------------------------
            ! Inertial force from UpdateStates (Includes a few other terms), root frame
            CASE (BldNd_FFixr,BldNd_FFiyr,BldNd_FFizr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fi(1:3,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_FFixr)   ! Gyroscopic force Fc x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_FFiyr)   ! Gyroscopic force Fc y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_FFizr)   ! Gyroscopic force Fc z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO
            CASE (BldNd_MFixr,BldNd_MFiyr,BldNd_MFizr)
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes      ! Index to current output
                  nelem = p%OutNd2NdElem(2,idx_node)
                  idx_node_in_elem = p%OutNd2NdElem(1,idx_node)
                  temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1), MATMUL(OtherState%GlbRot,m%qp%Fi(4:6,idx_node_in_elem,nelem)))
                  SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
                     CASE (BldNd_MFixr)   ! Gyroscopic moment Fc about x, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(1)
                     CASE (BldNd_MFiyr)   ! Gyroscopic moment Fc about y, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(2)
                     CASE (BldNd_MFizr)   ! Gyroscopic moment Fc about z, root frame
                        y%WriteOutput( IdxOutList ) = temp_vec(3)
                  END SELECT
               ENDDO



            CASE DEFAULT
               CALL SetErrStat( ErrID_Severe, "Coding error.  Output channel not properly set.",ErrStat,ErrMsg,RoutineName )
               DO idx_node=1,y%BldMotion%NNodes    ! Note p%node_total is total number of nodes including all elements
                  IdxOutList = p%NumOuts + idx_node + (IdxChan-1)*y%BldMotion%NNodes
                  y%WriteOutput( IdxOutList ) = 0.0_ReKi
               ENDDO

         END SELECT

   ENDDO ! Loop over the output channel list


END SUBROUTINE Calc_WriteBldNdOutput

!..................................................................................................................................
SUBROUTINE BldNdOuts_SetParameters(InitInp, InputFileData, p, ErrStat, ErrMsg)
   type(BD_InitInputType),       intent(in   )  :: InitInp           !< Input data for initialization routine
   type(BD_InputFile),           intent(in   )  :: InputFileData     !< data from the input file
   type(BD_ParameterType),       intent(inout)  :: p                 !< Parameters  ! intent(out) only because it changes p%NdIndx
   integer(IntKi),               intent(  out)  :: ErrStat           !< Error status of the operation
   character(*),                 intent(  out)  :: ErrMsg            !< Error message if ErrStat /= ErrID_None

   integer(IntKi)                               :: i
   integer(IntKi)                               :: ErrStat2
   character(ErrMsgLen)                         :: ErrMsg2
   CHARACTER(*), PARAMETER                      :: RoutineName = ('BldNdOuts_SetParameters')


      ! Initialize values
   ErrStat = ErrID_None
   ErrMsg = ""


      ! Set the parameter to store number of requested Blade Node output sets
   p%BldNd_NumOuts = InputFileData%BldNd_NumOuts

   IF ( p%BldNd_NumOuts == 0 ) THEN
   
      p%BldNd_TotNumOuts = 0  ! default to no nodal outputs
      
   ELSE
   
      ! Check if the blade node array to output is valid: p%BldNd_BlOutNd
      ! TODO: this value is not read in by the input file reading yet, so setting to all blade nodes
      !        -- check if list handed in is of nodes that exist (not sure this is ever checked)
      !        -- Make sure the nodes actually exist on the y%BldMotion mesh
      !        -- Sort the order of the list handed in
      !        -- copy values over


         ! Temporary workaround here:
      ALLOCATE ( p%BldNd_BlOutNd( size(p%NdIndxInverse) ) , STAT=ErrStat2 )
         IF ( ErrStat2 /= 0_IntKi )  THEN
            CALL SetErrStat( ErrID_Fatal,"Error allocating memory for the BeamDyn BldNd_BlOutNd array.", ErrStat, ErrMsg, RoutineName )
            RETURN
         ENDIF
      DO I=1,size(p%BldNd_BlOutNd)
         p%BldNd_BlOutNd(i) = i
      ENDDO

            ! Set the total number of outputs ( requested channel groups * number requested nodes )
      p%BldNd_TotNumOuts = p%BldNd_NumOuts * SIZE(p%BldNd_BlOutNd)
   
      call BldNdOuts_SetOutParam(InputFileData%BldNd_OutList, p, ErrStat2, ErrMsg2 )
         call setErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
         if (ErrStat >= AbortErrLev) return

   END IF


END SUBROUTINE BldNdOuts_SetParameters
!**********************************************************************************************************************************
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these 
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary. 
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine checks to see if any requested output channel names (stored in the OutList(:)) are invalid. It returns a 
!! warning if any of the channels are not available outputs from the module.
!!  It assigns the settings for OutParam(:) (i.e, the index, name, and units of the output channels, WriteOutput(:)).
!!  the sign is set to 0 if the channel is invalid.
!! It sets assumes the value p%NumOuts has been set before this routine has been called, and it sets the values of p%OutParam here.
!! 
!! This routine was generated by Write_ChckOutLst.m using the parameters listed in OutListParameters.xlsx at 12-Dec-2017 20:48:14.
SUBROUTINE BldNdOuts_SetOutParam(BldNd_OutList, p, ErrStat, ErrMsg )
!..................................................................................................................................

   IMPLICIT                        NONE

      ! Passed variables

   CHARACTER(ChanLen),        INTENT(IN)     :: BldNd_OutList(:)                  !< The list out user-requested outputs
   TYPE(BD_ParameterType),    INTENT(INOUT)  :: p                                 !< The module parameters
   INTEGER(IntKi),            INTENT(OUT)    :: ErrStat                           !< The error status code
   CHARACTER(*),              INTENT(OUT)    :: ErrMsg                            !< The error message, if an error occurred

      ! Local variables

   INTEGER                      :: ErrStat2                                        ! temporary (local) error status
   INTEGER                      :: I                                               ! Generic loop-counting index
   INTEGER                      :: J                                               ! Generic loop-counting index
   INTEGER                      :: INDX                                            ! Index for valid arrays

   LOGICAL                      :: InvalidOutput(1:BldNd_MaxOutPts)                      ! This array determines if the output channel is valid for this configuration
   CHARACTER(ChanLen)           :: OutListTmp                                      ! A string to temporarily hold OutList(I)
   CHARACTER(*), PARAMETER      :: RoutineName = "BldNdOuts_SetOutParam"

   CHARACTER(OutStrLenM1), PARAMETER  :: ValidParamAry(177) =  (/  &   ! This lists the names of the allowed parameters, which must be sorted alphabetically
                               "ABSXG   ","ABSXR   ","ABSYG   ","ABSYR   ","ABSZG   ", &
                               "ABSZR   ","AXB     ","AXL     ","AYB     ","AYL     ","AZB     ","AZL     ", &
                               "DFXL    ","DFXR    ","DFYL    ","DFYR    ","DFZL    ","DFZR    ","DMXL    ", &
                               "DMXR    ","DMYL    ","DMYR    ","DMZL    ","DMZR    ","FFBXL   ","FFBXR   ","FFBYL   ", &
                               "FFBYR   ","FFBZL   ","FFBZR   ","FFCXL   ","FFCXR   ","FFCYL   ","FFCYR   ","FFCZL   ", &
                               "FFCZR   ","FFDXL   ","FFDXR   ","FFDYL   ","FFDYR   ","FFDZL   ","FFDZR   ","FFGXL   ", &
                               "FFGXR   ","FFGYL   ","FFGYR   ","FFGZL   ","FFGZR   ","FFIXL   ","FFIXR   ","FFIYL   ", &
                               "FFIYR   ","FFIZL   ","FFIZR   ","FXB     ","FXL     ","FXR     ","FYB     ","FYL     ", &
                               "FYR     ","FZB     ","FZL     ","FZR     ","MFBXL   ","MFBXR   ","MFBYL   ","MFBYR   ", &
                               "MFBZL   ","MFBZR   ","MFCXL   ","MFCXR   ","MFCYL   ","MFCYR   ","MFCZL   ","MFCZR   ", &
                               "MFDXL   ","MFDXR   ","MFDYL   ","MFDYR   ","MFDZL   ","MFDZR   ","MFGXL   ","MFGXR   ", &
                               "MFGYL   ","MFGYR   ","MFGZL   ","MFGZR   ","MFIXL   ","MFIXR   ","MFIYL   ","MFIYR   ", &
                               "MFIZL   ","MFIZR   ","MXB     ","MXL     ","MXR     ","MYB     ","MYL     ","MYR     ", &
                               "MZB     ","MZL     ","MZR     ","PFXL    ","PFYL    ","PFZL    ","PMXL    ","PMYL    ", &
                               "PMZL    ","PXB     ","PXG     ","PYB     ","PYG     ","PZB     ","PZG     ","QXB     ", &
                               "QXL     ","QYB     ","QYL     ","QZB     ","QZL     ","RAXL    ","RAXR    ","RAYL    ", &
                               "RAYR    ","RAZL    ","RAZR    ","RDXR    ","RDYR    ", &
                               "RDZR    ","RVXG    ","RVXL    ","RVXR    ","RVYG    ","RVYL    ","RVYR    ","RVZG    ", &
                               "RVZL    ","RVZR    ","TAXL    ","TAXR    ","TAYL    ", &
                               "TAYR    ","TAZL    ","TAZR    ","TDXR    ","TDYR    ","TDZR    ", &
                               "TVXG    ","TVXL    ","TVXR    ","TVYG    ","TVYL    ","TVYR    ","TVZG    ","TVZL    ", &
                               "TVZR    ","UXB     ","UYB     ","UZB     ","VXB     ","VXG     ", &
                               "VXL     ","VYB     ","VYG     ","VYL     ","VZB     ","VZG     ","VZL     ","WXB     ", &
                               "WXG     ","WXL     ","WYB     ","WYG     ","WYL     ","WZB     ","WZG     ","WZL     "/)
   INTEGER(IntKi), PARAMETER :: ParamIndxAry(177) =  (/ &                            ! This lists the index into AllOuts(:) of the allowed parameters ValidParamAry(:)
                                 BldNd_AbsXg ,    BldNd_AbsXr ,    BldNd_AbsYg ,    BldNd_AbsYr ,    BldNd_AbsZg , &
                                   BldNd_AbsZr ,     BldNd_TAxr ,     BldNd_TAxl ,     BldNd_TAyr ,     BldNd_TAyl ,     BldNd_TAzr ,     BldNd_TAzl ,  &
                                  BldNd_DFxL ,     BldNd_DFxR ,     BldNd_DFyL ,     BldNd_DFyR ,     BldNd_DFzL ,     BldNd_DFzR ,     BldNd_DMxL , &
                                    BldNd_DMxR ,     BldNd_DMyL ,     BldNd_DMyR ,     BldNd_DMzL ,     BldNd_DMzR ,    BldNd_FFbxl ,    BldNd_FFbxr ,    BldNd_FFbyl , &
                                   BldNd_FFbyr ,    BldNd_FFbzl ,    BldNd_FFbzr ,    BldNd_FFcxl ,    BldNd_FFcxr ,    BldNd_FFcyl ,    BldNd_FFcyr ,    BldNd_FFczl , &
                                   BldNd_FFczr ,    BldNd_FFdxl ,    BldNd_FFdxr ,    BldNd_FFdyl ,    BldNd_FFdyr ,    BldNd_FFdzl ,    BldNd_FFdzr ,    BldNd_FFgxl , &
                                   BldNd_FFgxr ,    BldNd_FFgyl ,    BldNd_FFgyr ,    BldNd_FFgzl ,    BldNd_FFgzr ,    BldNd_FFixl ,    BldNd_FFixr ,    BldNd_FFiyl , &
                                   BldNd_FFiyr ,    BldNd_FFizl ,    BldNd_FFizr ,      BldNd_Fxr ,      BldNd_FxL ,      BldNd_Fxr ,      BldNd_Fyr ,      BldNd_FyL , &
                                     BldNd_Fyr ,      BldNd_Fzr ,      BldNd_FzL ,      BldNd_Fzr ,    BldNd_MFbxl ,    BldNd_MFbxr ,    BldNd_MFbyl ,    BldNd_MFbyr , &
                                   BldNd_MFbzl ,    BldNd_MFbzr ,    BldNd_MFcxl ,    BldNd_MFcxr ,    BldNd_MFcyl ,    BldNd_MFcyr ,    BldNd_MFczl ,    BldNd_MFczr , &
                                   BldNd_MFdxl ,    BldNd_MFdxr ,    BldNd_MFdyl ,    BldNd_MFdyr ,    BldNd_MFdzl ,    BldNd_MFdzr ,    BldNd_MFgxl ,    BldNd_MFgxr , &
                                   BldNd_MFgyl ,    BldNd_MFgyr ,    BldNd_MFgzl ,    BldNd_MFgzr ,    BldNd_MFixl ,    BldNd_MFixr ,    BldNd_MFiyl ,    BldNd_MFiyr , &
                                   BldNd_MFizl ,    BldNd_MFizr ,      BldNd_Mxr ,      BldNd_MxL ,      BldNd_Mxr ,      BldNd_Myr ,      BldNd_MyL ,      BldNd_Myr , &
                                     BldNd_Mzr ,      BldNd_MzL ,      BldNd_Mzr ,     BldNd_PFxL ,     BldNd_PFyL ,     BldNd_PFzL ,     BldNd_PMxL ,     BldNd_PMyL , &
                                    BldNd_PMzL ,    BldNd_AbsXr ,    BldNd_AbsXg ,    BldNd_AbsYr ,    BldNd_AbsYg ,    BldNd_AbsZr ,    BldNd_AbsZg ,     BldNd_RAxr , &
                                    BldNd_RAxl ,     BldNd_RAyr ,     BldNd_RAyl ,     BldNd_RAzr ,     BldNd_RAzl ,     BldNd_RAxl ,     BldNd_RAxr ,     BldNd_RAyl , &
                                    BldNd_RAyr ,     BldNd_RAzl ,     BldNd_RAzr ,     BldNd_RDxr ,     BldNd_RDyr , &
                                    BldNd_RDzr ,     BldNd_RVxg ,     BldNd_RVxl ,     BldNd_RVxr ,     BldNd_RVyg ,     BldNd_RVyl ,     BldNd_RVyr ,     BldNd_RVzg , &
                                    BldNd_RVzl ,     BldNd_RVzr ,     BldNd_TAxl ,     BldNd_TAxr ,     BldNd_TAyl , &
                                    BldNd_TAyr ,     BldNd_TAzl ,     BldNd_TAzr ,     BldNd_TDxr ,     BldNd_TDyr ,     BldNd_TDzr , &
                                    BldNd_TVxg ,     BldNd_TVxl ,     BldNd_TVxr ,     BldNd_TVyg ,     BldNd_TVyl ,     BldNd_TVyr ,     BldNd_TVzg ,     BldNd_TVzl , &
                                    BldNd_TVzr ,     BldNd_TDxr ,     BldNd_TDyr ,     BldNd_TDzr ,     BldNd_TVxr ,     BldNd_TVxg , &
                                    BldNd_TVxl ,     BldNd_TVyr ,     BldNd_TVyg ,     BldNd_TVyl ,     BldNd_TVzr ,     BldNd_TVzg ,     BldNd_TVzl ,     BldNd_RVxr , &
                                    BldNd_RVxg ,     BldNd_RVxl ,     BldNd_RVyr ,     BldNd_RVyg ,     BldNd_RVyl ,     BldNd_RVzr ,     BldNd_RVzg ,     BldNd_RVzl /)
   CHARACTER(ChanLen), PARAMETER :: ParamUnitsAry(177) =  (/  &  ! This lists the units corresponding to the allowed parameters
                               "(m)      ","(m)      ","(m)      ","(m)      ","(m)      ", &
                               "(m)      ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ",&
                               "(N/m)    ","(N/m)    ","(N/m)    ","(N/m)    ","(N/m)    ","(N/m)    ","(N-m/m)  ", &
                               "(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N)      ","(N)      ","(N)      ", &
                               "(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ", &
                               "(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ", &
                               "(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ", &
                               "(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ","(N)      ", &
                               "(N)      ","(N)      ","(N)      ","(N)      ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N)      ","(N)      ","(N)      ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(m)      ","(m)      ","(m)      ","(m)      ","(m)      ","(m)      ","(deg/s^2)", &
                               "(deg/s^2)","(deg/s^2)","(deg/s^2)","(deg/s^2)","(deg/s^2)","(deg/s^2)","(deg/s^2)","(deg/s^2)", &
                               "(deg/s^2)","(deg/s^2)","(deg/s^2)","(-)      ","(-)      ", &
                               "(-)      ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ", &
                               "(deg/s)  ","(deg/s)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ", &
                               "(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m)      ","(m)      ","(m)      ", &
                               "(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ", &
                               "(m/s)    ","(m)      ","(m)      ","(m)      ","(m/s)    ","(m/s)    ", &
                               "(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(m/s)    ","(deg/s)  ", &
                               "(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(deg/s)  "/)


      ! Initialize values
   ErrStat = ErrID_None
   ErrMsg = ""
   InvalidOutput = .FALSE.


!   ..... Developer must add checking for invalid inputs here: .....

      ! these outputs are not valid for FE (Gauss)
   IF (p%BldMotionNodeLoc==BD_MESH_FE) then
      InvalidOutput( BldNd_FFbxl ) = .true.
      InvalidOutput( BldNd_FFbyl ) = .true.
      InvalidOutput( BldNd_FFbzl ) = .true.
      InvalidOutput( BldNd_MFbxl ) = .true.
      InvalidOutput( BldNd_MFbyl ) = .true.
      InvalidOutput( BldNd_MFbzl ) = .true.
      InvalidOutput( BldNd_FFcxl ) = .true.
      InvalidOutput( BldNd_FFcyl ) = .true.
      InvalidOutput( BldNd_FFczl ) = .true.
      InvalidOutput( BldNd_MFcxl ) = .true.
      InvalidOutput( BldNd_MFcyl ) = .true.
      InvalidOutput( BldNd_MFczl ) = .true.
      InvalidOutput( BldNd_FFdxl ) = .true.
      InvalidOutput( BldNd_FFdyl ) = .true.
      InvalidOutput( BldNd_FFdzl ) = .true.
      InvalidOutput( BldNd_MFdxl ) = .true.
      InvalidOutput( BldNd_MFdyl ) = .true.
      InvalidOutput( BldNd_MFdzl ) = .true.
      InvalidOutput( BldNd_FFgxl ) = .true.
      InvalidOutput( BldNd_FFgyl ) = .true.
      InvalidOutput( BldNd_FFgzl ) = .true.
      InvalidOutput( BldNd_MFgxl ) = .true.
      InvalidOutput( BldNd_MFgyl ) = .true.
      InvalidOutput( BldNd_MFgzl ) = .true.
      InvalidOutput( BldNd_FFbxr ) = .true.
      InvalidOutput( BldNd_FFbyr ) = .true.
      InvalidOutput( BldNd_FFbzr ) = .true.
      InvalidOutput( BldNd_MFbxr ) = .true.
      InvalidOutput( BldNd_MFbyr ) = .true.
      InvalidOutput( BldNd_MFbzr ) = .true.
      InvalidOutput( BldNd_FFcxr ) = .true.
      InvalidOutput( BldNd_FFcyr ) = .true.
      InvalidOutput( BldNd_FFczr ) = .true.
      InvalidOutput( BldNd_MFcxr ) = .true.
      InvalidOutput( BldNd_MFcyr ) = .true.
      InvalidOutput( BldNd_MFczr ) = .true.
      InvalidOutput( BldNd_FFdxr ) = .true.
      InvalidOutput( BldNd_FFdyr ) = .true.
      InvalidOutput( BldNd_FFdzr ) = .true.
      InvalidOutput( BldNd_MFdxr ) = .true.
      InvalidOutput( BldNd_MFdyr ) = .true.
      InvalidOutput( BldNd_MFdzr ) = .true.
      InvalidOutput( BldNd_FFgxr ) = .true.
      InvalidOutput( BldNd_FFgyr ) = .true.
      InvalidOutput( BldNd_FFgzr ) = .true.
      InvalidOutput( BldNd_MFgxr ) = .true.
      InvalidOutput( BldNd_MFgyr ) = .true.
      InvalidOutput( BldNd_MFgzr ) = .true.
      InvalidOutput( BldNd_FFixl ) = .true.
      InvalidOutput( BldNd_FFiyl ) = .true.
      InvalidOutput( BldNd_FFizl ) = .true.
      InvalidOutput( BldNd_MFixl ) = .true.
      InvalidOutput( BldNd_MFiyl ) = .true.
      InvalidOutput( BldNd_MFizl ) = .true.
      InvalidOutput( BldNd_FFixr ) = .true.
      InvalidOutput( BldNd_FFiyr ) = .true.
      InvalidOutput( BldNd_FFizr ) = .true.
      InvalidOutput( BldNd_MFixr ) = .true.
      InvalidOutput( BldNd_MFiyr ) = .true.
      InvalidOutput( BldNd_MFizr ) = .true.
   END IF
   IF (.NOT. (p%OutInputs .and. p%BldMotionNodeLoc/=BD_MESH_FE)) then
      ! Distributed output channels not allowed yet as the mapping of the m%u_DistrLoad_at_y only exists for MESH_QP with p%OutInputs set to true
      InvalidOutput( BldNd_DFxL  ) = .true.
      InvalidOutput( BldNd_DFyL  ) = .true.
      InvalidOutput( BldNd_DFzL  ) = .true.
      InvalidOutput( BldNd_DMxL  ) = .true.
      InvalidOutput( BldNd_DMyL  ) = .true.
      InvalidOutput( BldNd_DMzL  ) = .true.
      InvalidOutput( BldNd_DFxr  ) = .true.
      InvalidOutput( BldNd_DFyr  ) = .true.
      InvalidOutput( BldNd_DFzr  ) = .true.
      InvalidOutput( BldNd_DMxr  ) = .true.
      InvalidOutput( BldNd_DMyr  ) = .true.
      InvalidOutput( BldNd_DMzr  ) = .true.
   END IF
   
!   ................. End of validity checking .................


   !-------------------------------------------------------------------------------------------------
   ! Allocate and set index, name, and units for the output channels
   ! If a selected output channel is not available in this module, set error flag.
   !-------------------------------------------------------------------------------------------------

   ALLOCATE ( p%BldNd_OutParam(1:p%BldNd_NumOuts) , STAT=ErrStat2 )
   IF ( ErrStat2 /= 0_IntKi )  THEN
      CALL SetErrStat( ErrID_Fatal,"Error allocating memory for the BeamDyn BldNd_OutParam array.", ErrStat, ErrMsg, RoutineName )
      RETURN
   ENDIF


      ! Set index, name, and units for all of the output channels.
      ! If a selected output channel is not available by this module set ErrStat = ErrID_Warn.

   DO I = 1,p%BldNd_NumOuts

      p%BldNd_OutParam(I)%Name  = "_"//BldNd_OutList(I)
      OutListTmp          = BldNd_OutList(I)
      p%BldNd_OutParam(I)%SignM = 1   ! this won't be used

      CALL Conv2UC( OutListTmp )    ! Convert OutListTmp to upper case


      Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry )

      IF ( Indx > 0 ) THEN ! we found the channel name
         IF ( InvalidOutput( ParamIndxAry(Indx) ) ) THEN  ! but, it isn't valid for these settings
            p%BldNd_OutParam(I)%Indx  = 0                 ! pick any valid channel (I just picked "Time=0" here because it's universal)
            p%BldNd_OutParam(I)%Units = "INVALID"
            p%BldNd_OutParam(I)%SignM = 0
         ELSE
            p%BldNd_OutParam(I)%Indx  = ParamIndxAry(Indx)
            p%BldNd_OutParam(I)%Units = ParamUnitsAry(Indx) ! it's a valid output
         END IF
      ELSE ! this channel isn't valid
         p%BldNd_OutParam(I)%Indx  = 0                    ! pick any valid channel (I just picked "Time=0" here because it's universal)
         p%BldNd_OutParam(I)%Units = "INVALID"
         p%BldNd_OutParam(I)%SignM = 0                    ! multiply all results by zero

         CALL SetErrStat(ErrID_Fatal, TRIM(p%BldNd_OutParam(I)%Name)//" is not an available output channel.",ErrStat,ErrMsg,RoutineName)
      END IF

   END DO

   RETURN
END SUBROUTINE BldNdOuts_SetOutParam
!----------------------------------------------------------------------------------------------------------------------------------
!End of code generated by Matlab script
!**********************************************************************************************************************************

END MODULE BeamDyn_BldNdOuts_IO
